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ABSTRACT 

"<N"' 

o : 

Far Ultraviolet (FUV, 6 eV< hv <13.6 eV) radiation has been recognized as 
the main source of heating of the neutral interstellar gas, and, as a consequence, 
'-^ | it determines whether the thermal balance of the neutral gas results in cold 

[JLh ; (T ~ 50 - 100K) clouds (CNM), warm (T ~ W^K) clouds (WNM), or a combination 

of the two. High FUV fields convert the neutral gas to WNM, while low fields result 
in CNM. The knowledge of how these fractions depend on the FUV sources (i.e. the 
y—i . star formation rate, the IMF, and the size distribution of associations) is a basic step 

in building any detailed model of the large scale behavior of the ISM and the mutual 
| relation between the ISM and the star formation rate in a galaxy. 

The sources of FUV radiation are the short-lived massive stars that generally 
originate in associations that form in Giant Molecular Clouds present in the galactic 
disk. Using McKee & Williams' (1997) distribution of birthrates for OB associations 
■ in the Galaxy, we determine the expected behavior of the time-dependent FUV 

field for random positions in the local ISM. The FUV field is calculated in two 
bands (912 - 1100 A and 912 - 2070 A) and at the wavelength 1400 A. In terms of 
U-i7 = U /(10 -17 erg cm' 3 A -1 ), where U is the energy density of the radiation field in 
some band, we find (mean, median) values at the solar circle of U—n =(15.7, 7.4) and 
(14.2, 7.2) for the [912-1100 A] and [912-2070 A] bands, respectively. At 1400 A we find 
J^j ■ (mean, median) values of U-n =(14.4, 7.5). Our median value for the [912-2070 A] 

band is Go = 1.6 times Habing's (1968) value for the radiation field at the solar circle 
in this band, and quite close to Draine's (1976) value, Go = 1.7. Both the latter values 
are based on observations of sources of FUV radiation in the solar neighborhood, so all 
three values are close to observed values. Due to attenuation by dust, only associations 
within about 500 pc contribute significantly to the energy density at a given point. 
Large angle scattering produces a diffuse field that is about 10% of the field produced 
by the sum of direct and small angle (< 5°) scattering from discrete sources (the 
associations), as observed. At a point exposed to the median radiation field, the 
brightest association typically produces about 20% of the total energy density. At a 
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point exposed to an above average radiation field, the brightest association produces 
most of the energy density. Therefore, the FUV field is asymmetric at a given point, 
and the asymmetry grows for higher fields. 

The FUV field fluctuates with a variety of amplitudes, the larger ones being less 
frequent. The mean field is about twice the median field because of these fluctuations, 
or spikes, in the radiation field. These spikes, which last ~ 30 Myr, are caused by the 
infrequent birth of nearby associations. For spikes that are significantly higher than 
the mean field, the time interval between spikes is ~ 2£7_ 15 Gyr. We also model shorter 
duration spikes caused by runaway OB stars. The presence of a fluctuating heating 
rate created by the fluctuating FUV field converts CNM to WNM and vice versa. 

Subject headings: ISM: structure — ISM: evolution — Stars: formation — Stars: mass 
function 

1. Introduction 

Neutral interstellar gas is the dominant mass component of the interstellar medium (ISM), 
and the atomic gas exists as two phases in rough pressure equilibrium: the warm neutral medium 
(WNM) with hydrogen densities at the solar circle of n ~ 0.3 cm -3 and T ~ 8000 K, and the 
cold neutral medium (CNM) with n ~ 40 cm" 3 and T ~ 70 K (Wolfire et al. 1995, 2002). With 
sufficient shielding column density, N > 10 20 ~ 21 cm -2 , the CNM clouds harbor molecular interiors, 
and for N > 10 22 cm -2 they become gravitationally bound and stars may form in their molecular 
interiors (e.g., McKee 1989). Most of the star formation in the Galaxy occurs in Giant Molecular 
Clouds (GMCs) with large masses M ~ 10 5-6 M Q and columns N ~ 10 22 cm~ 2 . 

In spiral galaxies FUV radiation from OB stars heats, photodissociates, and controls the state 
of the neutral interstellar gas. High FUV fields heat CNM to WNM; low FUV fields allow WNM to 
cool to CNM (Wolfire et al. 1995, 2002). If, as seems likely, the formation rate of stars in galaxies 
depends on the mass fraction of CNM gas, then a possible feedback process occurs between stars 
and the ISM that may control the global rate of star formation in a galaxy (Parravano 1988, 1989, 
1991; Corbelli & Salpeter 1995). High star formation rates lead to large numbers of OB stars that 
create high FUV fluxes that destroy the sources of star formation. Low rates of star formation, 
on the other hand, may lead to enhanced growth rates of star-forming clouds and thereby to 
enhanced star formation rates. In this paper we focus on the calculation of the fluctuating FUV 
field in the local region of the Milky Way Galaxy, deferring to future papers the study of the 
fractional amounts of CNM and WNM produced by the FUV field, and the possible regulation of 
star formation rates in galaxies. 

OB stars are generally born in associations that form inside GMCs present in the galactic 
disk. The FUV field at a given point in a galaxy then depends largely on the distribution of 
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stellar masses in relatively nearby associations and on the birth-rate and size distribution of these 
associations. The FUV luminosity of an association is an increasing function of its size, but the 
birth-rate of associations decreases with size. This fact, coupled with the fact that the probability 
of a young association existing within a distance r of a point decreases rapidly with decreasing r, 
leads to a fluctuating FUV field, with the larger fluctuations being less frequent. The abundant 
low luminosity or distant associations help to build an FUV background, over which is added 
bursts of FUV flux coming from relatively rare events (the birth of nearby large associations) that 
can significantly increase the FUV field for about ~ 30 to 50 Myr. This period corresponds to the 
time that an association produces stars (~ 20 Myr, generally in a mode of several bursts) plus the 
time that early B stars, which dominate the FUV production, remain alive after the end of the 
star formation process in the GMC. The time variation of the FUV field at any given point in a 
galaxy maps directly onto the spatial variation of the FUV field in the galaxy at any instant of 
time: those regions very close to luminous associations will have exceptionally high FUV fluxes, 
whereas large volumes of space will have much lower FUV fluxes. 

The presence of a time and space-dependent FUV field has a major impact on the structure 
of the interstellar medium. In a subsequent paper we calculate the response of interstellar gas 
to this fluctuating FUV field, estimate the average fractions of the gas in the CNM and WNM 
phases, and determine the typical time that the gas remains in each phase. The knowledge of these 
quantities and how they depend on the gas properties (i.e., mean density and composition) and on 
the FUV-sources (i.e., star formation rate SFR, initial mass function IMF, and size distribution of 
associations) is a basic step in building any detailed model of the large scale behavior of the ISM 
and the mutual relation between the ISM and the SFR. 

FUV photons heat the neutral interstellar gas primarily by photoejecting energetic (~ 1 
eV) electrons into the gas from interstellar dust grains (Watson, 1972; Jura, 1976; de Jong, 
1977; Draine, 1978) and large stable molecules like polycyclic aromatic hydrocarbons (PAHs) 
(d'Hendecourt and Leger, 1987; Lepp and Dalgarno, 1988; Verstraete et al., 1990; Bakes and 
Tielens, 1994). The relevant wavelength range for this "grain photoelectric heating" process is 
determined by the frequency dependence of the absorption cross sections and the photoelectric 
yields of grains and PAHs. Since we are interested in the neutral ISM, the Lyman threshold 
(912 A) is the shortest wavelength of the relevant band. The long wavelength limit is somewhat 
arbitrary. It depends on the wavelength at which both the absorption cross section and the yield 
has fallen significantly for grains as well as for PAHs (Draine 1978; Bakes & Tielens 1994, and 
references therein). In general, the band between 6 eV< hv < 13.6 eV (i.e. [912 — 2070 A]) has 
been considered as the primary band for grain photoelectric heating, and a mean value of the flux 
in the band is often used as a parameter to calculate the thermochemical equilibrium of the neutral 
gas. However, for a given total flux in the band, a change in its energy distribution can have a 
large effect on the photodissociation rates and even the heating rates. In particular, the main 
processes of H2 and CO destruction is photodissociation by radiation in the band [912 — 1100 A]. 
We therefore follow two primary bands: 
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a) The band [912 - 2070 A], hereafter the "FUV band". The mean energy density in this 
band is related to the parameter Go that is often used to parametrize the photoelectric dust 
heating rate (Tielens & Hollenbach 1985). Let {/band = {/(band) / Ab an d be the average radiation 
energy density per A in the wavelength band Aband- The energy density of FUV radiation is 
generally expressed in units of the Habing (1968) typical energy density at the solar circle averaged 
between 6 eV< hv < 13.6 eV (i.e. [912 - 2070 A]), 

„ _ {/fuv {/fuv 

° ~ ^FUV ~ 46 x 10-17 erg cm" 3 A" 1 

If the radiation field is one-dimensional, it is convenient to express the Habing field as a flux, 
c{/^ uv A F uv = 1-6 x 10~ 3 ergcm^s -1 . 

b) The band [912 — 1100 A], hereafter the " H2 band". The flux in this band determines 
the dissociation rate of H2 and CO. It is also useful to consider separately the contribution of this 
band to the photoelectic heating. Even though the H2 band covers only about 10% of the whole 
FUV wavelength band, both the absorption cross section and the yield of dust are significantly 
higher in the H2 band, so that it contributes disproportionately to the gas heating. On the other 
hand, the extinction in the H2 band is about twice the extinction in the FUV band, which reduces 
the importance of the H2 band heating. 

In addition to these two rather broad bands, simulations in a 1 A band around 1400 A are 
performed in order to compare the averages of the time dependent FUV model with observations. 
Many observations have been made of the interstellar field at 1400 A, although the observations 
are averages over bands considerably broader than 1 A. 

In our model, we include both direct and scattered FUV, and therefore, at a given time 
the local value of the FUV interstellar radiation field (hereafter FUV-ISRF) can be reasonably 
determined from the location of the FUV sources and the distribution and properties of the 
absorbing and scattering dust particles. In order to construct a time dependent model for FUV 
interstellar radiation we adopt: a) an initial mass function for stars in associations (§2.1); b) a 
stellar-mass dependence of the FUV and H2 band luminosities, and the main sequence lifetime 
for stars of various masses (§2.2); c) the dependence of the birth-rate of associations on the 
number of stars in the association, the history of the star formation rate in an association, and 
the distribution in the galactic plane of the association birthplaces from McKee & Williams 
(1997; hereafter MW97) (§3); and d) a vertical distribution of associations and dust (§4). We 
also present in §4 model simulations of the time dependence of the FUV-ISRF, together with the 
corresponding asymptotic values. The typical FUV-ISRF is computed as the median value in the 
time-sequence. We compare in §5 the mean and median values derived with the time-dependent 
model with previous estimations of the local field. In §6 we quantify the effect of runaway OB 
stars on the detailed time evolution of the FUV-ISRF. Finally, we summarize the results in §7. 
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2. Stellar Properties 

2.1. Stellar IMF and PDMF 

Let A/"*(> m) be the number of stars with masses greater than mM Q , and let 
dj\f*{m) = — d/V*(> m) be the number of stars with masses between m and m + dm. We 
express the present day mass function (PDMF) as dN*{m) / d\n m, which is simply the number of 
stars in a logarithmic mass interval. We denote the PDMF per unit area of Galactic disk by 

"* {m) = Z4dhW (2) 

This quantity is related to the quantity (p ms defined by Miller h Scalo (1979; hereafter MS79) by 
<T*(m) = ms log e = 0.434^ ms . 

The stellar birthrate as a function of mass — the IMF — is given by dN*(rn)/dlnm; note that 
A/"*(m) includes only stellar birth, not stellar death. The IMF per unit area of Galactic disk is 
<j*(m), which is related to MS79's £(logm) by / q(m,t)dt = 0.434£(log m). Sometimes the IMF 
is described by the fraction of stars above mass m [e.g., F n (> m) — MS79], but this depends on 
the uncertain and possibly variable number of low- mass stars in the IMF. Since we are interested 
in the stars that produce the FUV radiation, we shall generally characterize the IMF in terms of 
quantities related to the high-mass end of the IMF. 

For massive stars, the IMF can be approximated as a cut-off power law that extends up to a 
mass m u (MW97), 

- Al*u (m< m u ), (3) 



dlnm \ m 

where J\f m is approximately equal to the birthrate of stars with masses between 0.5m u and m u (for 
T = 1 this relation is exact). Young star clusters have PDMFs that are of this form also, although 
the parameters m u and V that enter this distribution may differ from those in the IMF due to 
stellar mass loss and, eventually, to supernovae. The birthrate of stars more massive than m is 

A/*(> m) = — 



m 



1 



(4) 



In most associations, there are very few stars with masses near m u , so it is convenient to express 
the IMF in terms of N^, the birthrate of high-mass stars (those with masses m > m^), 

*»). J f,fD(^ r , (5) 



dlnm \(f>hJ\ m 

where 4>h = 1 — { r nh/ m u) T is a numerical parameter that is of order unity; indeed, for the fiducial 
values of the parameters that we adopt below, <ph = 0.974. 



The value of rrih is somewhat arbitrary; we adopt rrih = 8 since most ionizing photons 
are produced by stars above this mass and since this is the minimum mass for core-collapse 
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supernovae. The value of T is also somewhat uncertain. Van Buren (1983) carried out a careful 
study of the high stellar mass IMF of the Galaxy, and concluded that T = 1.03. A fit to the 
high-mass end of Scalo's (1986) IMF gives T ~ 1.5. In a study of OB associations in the northern 
Milky Way, Massey, Johnson, & DeGioia-Eastwood (1995) found a mean value T = 1.1 ± 0.1. In 
the star cluster R136 located in the LMC, Massey and Hunter (1998) found T ~ 1.35. While it 
remains controversial whether the IMF is a universal function (Scalo 1998), there is no evidence 
that the variations seen in the IMF are anything but statistical in nature. We adopt a single 
power-law IMF with T = 1.35 for stars more massive than 8 M@. As described below and in 
Appendix A, the probabilistic distribution used here allows for fluctuations around this power law 
in a given association. 

The existence of an upper limit to the mass function is also controversial; for example, Massey 
and Hunter (1998) suggest that the apparent cutoff is due to the sparcity of very massive stars. 
This issue was addressed in a general way by MW97. If the distribution of a quantity x can be 
approximated by a cutoff-power law of the form in equation (||) , 



then the expected number of objects with x > x u in the absence of a physical cutoff is M u /T. If 
this number is large, then the cutoff must be physically significant, whereas if it is ^ 1 then it 
is not. In R136, Massey and Hunter (1998) find 29 stars more massive than 75 -M©, but none 
more massive than 155M (based on the stellar data of Vacca, Garmany, & Shull 1996). Setting 
m u = 155 in this case, we find M u /T = 17 .4 3> 1. Hence, we conclude that the upper mass limit 
in R136 is physically significant. In the Galaxy, it is not clear that m u is as large as 155. In their 
survey of OB associations in the northern Milky Way, Massey et al. (1995) found only one star 
more massive than 120 M©, whereas many more would have been expected if the IMF continued 
on to higher masses. We shall adopt m u = 120 here. 



As they evolve, stars change their radii and effective temperatures and therefore their FUV 
luminosities. However, we are interested in the FUV radiation field from a large population of 
stars of different ages and in different associations, and therefore the use of a mean luminosity 
during the main sequence (MS) lifetime, i ms , is an adequate approximation. During the large but 
infrequent fluctuations of the FUV field that occur when a single nearby association dominates 
the radiation field, this approximation is not as good as in the quiescent period when associations 
of different ages contribute to the FUV field. However, even in a single association, several 
generations of stars can contribute to the total luminosity of the association, and the use of 
time-averaged FUV luminosities provides an adequate description of the total FUV luminosity of 
the association. Therefore, in our model the time variation of the total FUV luminosity of an 




(6) 



2.2. 



Stellar Sources of FUV Radiation 
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association is due to the birth and death of its massive stars, and not to the luminosity changes in 
individual stars during their MS lifetime. 

To compute the mean main-sequence luminosities we have used the set of stellar evolutionary 
tracks generally referred as the Padova tracks (Bruzual 1999, and references therein), with the 
stellar flux in the H2 and FUV bands kindly provided by G. Bruzual and J. Mateu (based on the 
library of synthetic stellar spectra compiled by Lejeune et al. 1997ab). Hereafter, we denote as 
Lfuv and Lh2 the mean main-sequence luminosities in the FUV band, [912 — 2070 A], and in the 
H2 band, [912 — 1100 A], respectively. For stars with solar metallicity, Figures la and lb show 
respectively the resulting mean luminosities in the band [1100 — 2070 A] (i.e. Lfuv — £#2) and in 
the H2 band (i.e. Lwi) as a function of the zero age main sequence stellar mass m in units of Mq. 
In Figures la and lb the triangles correspond to the calculated mean MS luminosities using the 
stellar evolutionary tracks, and the lines are the eight-segment power law used to approximate the 
mass dependence of the luminosity in the simulations. Table 1 shows the adopted power laws for 
Lfuv — L#2> Ljj2, S49 (the ionizing photon luminosity, see below), and the main sequence lifetime 
^ms(^)- The power law approximations for the luminosities are good to within 10%, and those for 
the main sequence lifetime are good to within 5%. 

When a massive star is young and embedded or near its natal cloud, the natal cloud absorbs 
a fraction of the emitted radiation. Due to their shorter lifetimes, more massive stars are more 
affected by this obscuration. We assume that the fraction of radiation locally absorbed goes 
linearly from 3/4 to zero in a time t bsc- in the nearest region of massive star formation, the Orion 
Molecular Cloud, the Trapezium generation of massive stars and the other Blaauw subgroups of 
massive stars apparently broke out of the obscuring cloud material in few million years. A more 
exact analysis can be done by quantitatively considering the breakout of OB stars from their 
clouds, as was done in the models of Whitworth (1979) and Williams h McKee (1997). The 
latter models indicate i bsc — 4 Myr, the number that we adopt. Because B stars are dominant 
contributors to the interstellar FUV field (see below), and because they live considerably longer 
than tobsc, the exact value of i b sc does not significantly affect our results, i.e., most of the FUV 
radiation is not absorbed locally in the natal clouds of the OB stars. 

In order to show the time averaged contribution of the different mass ranges to the FUV 
radiation field, Figures lc and Id show the energy output in each band for a mass distribution of 
stars following a Present Day Mass Function (PDMF) of MS stars. The PDMF <f*(m) has been 
constructed assuming a constant star formation rate (appropriate for the solar circle-see §2.4 
below), a single power-law for the IMF (§2.1), and weighting the stars by their main-sequence 
lifetimes. Rocha-Pinto et al. (2000) have presented evidence that the Galactic star formation 
rate changes on time scales of order 400 Myr, which corresponds to the main-sequence lifetime 
of a 2.9 M & star. Almost all the FUV radiation is produced by stars more massive than this, so 
such a time variation in the Galactic star formation rate would not affect our results. We note 
that the results in Figures lc and Id are useful only to make time-averaged estimates. To follow 
the evolution of the FUV-ISRF when the SFR is time and position dependent, it is necessary 
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to consider that obscuration acts during the period t bso starting at the moment when each 
generation of stars are formed. 

In the FUV band the corresponding mean emissivity per unit disk area is ~ 13ALq/pc 2 (i.e., 
the area below the heavy line curve in Fig. lc). After correcting for obscuration, the emissivity 
is ~ lO.9L /pc 2 (i.e. the area below the thin line curve in Fig. lc). Therefore, ~ 2.5-Lq/pc 2 
is absorbed by the natal clouds, representing ~ 19% of the total mean emissivity per unit disk 
area.0 The distribution of emissivity with stellar mass has a maximum at ~ 6M Q (i.e., B5), but 
only ~ 16% of the total emission is produced by stars with masses below this maximum. Half of 
the FUV radiation is produced by stars with m > 18, but half of the escaping FUV radiation is 
produced by stars with m > 13 (i.e., earlier than Bl). Stars with masses less than 3Mq only 
contribute ~ 3%. Therefore, the SFR history in the past ~ 350 Myr [i.e. i ms (3M )] is sufficient 
to simulate the present FUV radiation field. 

The characteristics of the FUV-emitting stars can be described in terms of their mean lifetime 
and luminosity. Ignoring absorption by the natal cloud, we calculate the luminosity- weighted 
lifetime as 

(t FUV ) = /^Fuv(ra)£ ms d\/; 

/ LFijy(m)dAf* 

which is 7.8 Myr for the stellar parameters we have adopted. By comparison, the corresponding 
lifetime for ionizing radiation is 4.0 Myr (§2.3). The mean FUV luminosity per high-mass star in 
the IMF is 



(LFw)h = "rr [ L FVY (m)dJ\f* = ( L FUV (m) (— ] dhxm. 
The mean FUV emissivity per unit disk area is related to the high-mass star formation rate as 



(8) 



Sfuv = ?*/i(-^FUv)h(*FUv)- (9) 

For the stellar parameters we have adopted (-Lfuv)/i = 5.8 x 10 4 L Q , and using the high-mass 
star formation rate derived in §2.4, Efuv = 13.4L /pc 2 in accordance with the value derived 
from the PDMF above. Since Spuv = 10.9Lq/pc 2 after correcting for obscuration, the 
obscuration-corrected mean luminosity is (Lfuv)/i = 4.7 x 10 L & . 

In the H2 band, the radiation field is dominated by more massive stars. The mean emissivity 
per unit disk area is ~ 3.2Lq/pc 2 . After correcting for obscuration, the emissivity is ~ 2.5L Q /pc 2 . 
The maximum contribution per unit mass occurs at m = 9 (i.e., B2 stars). After correcting 
for obscuration, about 85% of the total emission in the band is produced by stars with masses 



4 This calculation has interesting implications. If the diffuse ISM is optically thick to FUV photons, then 81% 
of the IR luminosity from galactic dust heated by FUV photons will arise from the diffuse gas. Therefore, when 
comparing IR observations of galaxies made with rather large beams, encompassing many GMCs and many FUV 
optical depths in the galactic plane, with PDR models (e.g., Malhotra et al. 2001, Kaufman et al. 1999) one will 
need to separate the IR continuum and [CII] and [OI] fine structure line emission into two components: one from the 
natal GMCs and one from the diffuse ISM. 
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m > 9, and stars with masses less than 5 M Q only contribute ~ 2% to the total emission (lifetime 
~ 100 Myr). Half of the total escaping emission is produced by stars with masses over ~ 21 Mq 
(lifetime of ~ 9 Myr) . Since the extinction in the H2 band is larger than the extinction in the FUV 
band, the radiation field in the H2 band is more dominated by nearby associations. Therefore, 
the amplitudes of the fluctuations in the H2 band are expected to be larger than in the FUV 
band, where more sources contribute to the field. For the stellar parameters we have adopted 
and ignoring absorption by the natal cloud, the luminosity- weighted lifetime is (£h 2 ) = 6.3 Myr 
and the mean H2 luminosity per high-mass star in the IMF is (Lh 2 )^ = 1.8 x 10 4 L Q . Since 
Eh 2 = 2.5L /pc 2 after correcting for obscuration, the obscuration-corrected mean luminosity is 
(L F iJv}h = 1-0 x 1O 4 L . 

Most of the conclusions in the preceding two paragraphs apply only to a distribution of FUV 
sources that is approximately constant in both space and time. In fact, the FUV sources are 
located in stellar associations that are spatially and temporally very inhomogeneous. As a result, 
the local SFR, and therefore the FUV and H2 band emissivities, are below their mean values most 
of the time (see §4.3). In addition, massive stars are scarce and therefore distant, so that they 
are subject to relatively large interstellar extinction. All of these effects will be included in our 
time dependent model. In any case, the results in Figures lc and Id show that it is sufficient to 
consider only the contribution of stars with masses larger than 2 or 3M Q to compute the time 
dependence of the FUV radiation field. 



2.3. Stellar Sources of Ionizing Radiation 

Ionizing radiation from hot stars is considered in the present model because we are calculating 
the FUV-ISRF in the neutral gas and therefore it is necessary to verify that the sampling point is 
outside of HII regions. In addition, we use the observations of the total ionizing photon production 
rate in the Galaxy to constrain the birth rate of associations (MW97; see §4) . 

Let s{m) be the time-averaged main-sequence ionizing photon luminosity (in photons s _1 ) of 
a star of initial mass m. Values of s{m) were kindly provided by Vacca (private communication) 
and are plotted in Figure 2a as S49 = s/(10 49 photons s -1 ). An analytic fit to these data is given 
in Table 1. In Vacca, Garmany, & Shull (1996) the ionizing luminosities are given as function of 
the spectroscopic mass and the evolutionary mass, but here the zero age main sequence mass is 
used because the IMF (§2.1) refers to this initial mass. 

The contribution of the different mass ranges to the ionizing photon production is shown in 
Figure 2b. Again a mass distribution of stars following a PDMF was assumed. For the adopted 
distribution (constant SFR and an IMF with power law V = 1.35, see §2.1), stars of about 30 M 
make the maximum contribution. Half of the ionizing photons are produced by stars above 
~ 55 M Q (lifetime ~ 3.8 Myr), and ~ 95% are produced by stars above ~ 20 M . 

The characteristics of the ionizing luminosity from a group of stars following an IMF that is a 
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cutoff power law can be described by two quantities (McKee 1989). The first is the mean ionizing 
photon luminosity per high-mass star, 

(s) h = 4r I s ( m )dK = T" / sH (—) dlnrn, (10) 

where the second equation follows from equation (||). The second quantity is the mean lifetime of 
an ionizing star, 

_ Js(m)t ms dN* 

\ r ion/ — - . s •- • (i-i-J 

J s(m)a/v* 

If each star produces Q(m) = s(m)t ms ionizing photons during its lifetime, then these equations 
imply that in a steady state the expected value of the ionizing luminosity from a group of stars is 



S = J Q{m)dN* = K h {s) h {t ion ). (12) 
For the stellar parameters we have adopted, (3)^49 = 0.42 and {t- lon ) = 4.0 x 10 6 yr, so that 

N*h = 6.0 x lO- 7 5 ^49, (13) 

where AC/i is measured in yr -1 and <ps allows for deviations from our adopted values of (s)h 
and (iion)- If fn u differs from our adopted value of 120, 4>s changes. For example, 4>s = 1.32 
for m u = 80 and <ps = 0.88 for m u = 150; a power-law approximation over this range gives 
(f>$ — (120/m M ) ' 65 . By comparison MW97, who adopted V = 1.5 and used a different set of stellar 
lifetimes but the same value of m u , found (s)h,i9 = 0.34 and (ii on ) = 3.7 x 10 6 yr, corresponding 
to <f> s = 1.32. 



2.4. Star Formation in the Galaxy 

What is the relation between the star formation rate locally and that in the Galaxy as a 
whole? Let <i*h(R) be the birthrate of massive stars (m > rrih) per unit area at Galactocentric 
radius R. MW97 found that the spatial distribution of massive star formation in the Galaxy 
(omitting the Galactic Center) can be approximated as s*h(R) on exp(—R/Hji) for 11 kpc > R > 
3 kpc, with Hr = 3.5 kpc. The total rate of massive star formation in the Galaxy is then 

/•ll kpc 

Kh,T = t*h 27reM(Ro-R)/H R ]R 2 dR (14) 

J 3 kpc 

= ^ h A e s, (15) 

where Rq = 8.5 kpc is the distance to the Galactic Center and <j#/j without an argument is the 
value of the massive star formation rate at the solar circle. [This definition of A e ff differs by a 
factor exp(i?o/^i?) from that in MW97]. Numerically, the effective area for massive star formation 
in the Galaxy is A c g = 530 kpc 2 . It should be borne in mind that star formation is observed to 
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occur beyond 11 kpc, but there are no giant radio HII regions, and therefore few if any large 
associations, there. 

Radio surveys of the Galaxy provide an estimate for the total ionizing photon production 
rate, St = 2.6 x 10 53 photons s _1 , good to a factor of about 1.5 (MW97). We use St to set 6^, 
since equations ( |l3| ) and ( |l5| ) yield 

U = 6.0 x l(T 7 0s f-j^) -» 29.405 kpc" 2 Myr" 1 , (16) 

where we have expressed the result in convenient units. This result is compared with other 
estimates in the literature in Figure 3. We note that the value <j*/i = 38 kpc -2 Myr -1 derived in 
MW97 was obtained by assuming T = 1.5 and adopting a different set of stellar lifetimes. 

In the Galaxy as a whole, the birthrate of massive stars is 6.0 x W^^sS^q^t = 0.015605 yr -1 . 
Correspondingly, the time between core-collapse supernovae in the Galaxy is tsN ~ (64/05) yr. 

The resulting PDMF for m > agrees reasonably well with the one derived by MS79 from 
observed data. Our PDMF is systematically higher than the one derived by Scalo (1986), although 
it is close to the top of his error bars. This trend can be explained by the fact that our PDMF 
corresponds to the mean star formation rate at the solar circle as determined from the Galactic 
distribution of emission of ionizing photons, whereas the observed local PDMF likely corresponds 
to the median star formation rate (see §4.3). Also, massive stars embedded or close to molecular 
clouds might be missing in the observed PDMF. 

3. Associations 

The sources of FUV radiation are the short-lived massive stars that generally form in 
associations. Therefore, these associations can be considered as point sources with an FUV 
luminosity that decreases with time as their massive members die. These associations tend to 
expand (typical velocities are ~ 1 — 2 km s^ 1 ), but this effect can be neglected because the 
associations remain small during the time that most of the FUV radiation is produced — recall 
that a large fraction of the FUV radiation that escapes from the natal cloud escapes during the 
first 16 Myr (§2.2). The fact that the main contributors to the FUV field are grouped has a 
tremendous effect on the temporal behavior of the radiation field. For example, the occurrence 
of a close association can raise the local radiation field by orders of magnitude until its massive 
star population dies. Therefore, we must model the birth-rate and size-distribution of associations 
in order to determine the time and space dependence of the interstellar FUV radiation field. 
We adopt the MW97 model, which explicitly gives the birthrate of associations as a function of 
the number of stars formed in the association. The birth-rate distribution in MW97 is derived 
from observations of HII regions, and reproduces the total ionizing luminosity, the local ionizing 
luminosity, and the rate of core collapse supernovae observed in the Galaxy. Moreover, the 
observed OB associations in the solar neighborhood (i.e., those that are young and large enough 
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to be identified) are consistent with the MW97 model, as we discuss below. The MW97 model is 
based on observations of the effects produced by Lyman continuum photons radiated by a small 
fraction of the stellar population (i.e. the very hot stars). A larger fraction of the star population 
contributes to the FUV radiation field; therefore an additional test to the MW97 model will be 
provided here. 



3.1. Evolution of the Luminosity of Associations 

The typical period of time that an association forms stars is of the order of 15 — 20Myr 
(Blaauw 1964, 1991; Heiles 1990). The rate at which stars form in an association is not continuous 
but occurs in bursts during the 15-20 Myr interval. We idealize this situation by assuming that 
each association produces five bursts, or generations, each of which has the same number of stars 
with a stellar mass distribution given by equation (||). We characterize the size of the association 
by A*/i, the number of massive stars that form over the 5 generations. The stars in each generation 
are assumed to form simultaneously. As in MW97, we assume that these five generations occur 
with a period of t g = 3.7 Myr. At any given time after the birth of an association, the luminosities 
£fuv ana - ^H 2 °f the entire association are the sum of the luminosities of the stars still living. 
However, as discussed in §2.2, not all of this UV luminosity can escape the parent cloud of the 
association. We estimated in §2.2 that for each generation the fraction of the radiation that 
escapes, / esc , goes linearly from 1/4 to one in a time t bsc ~ 4 Myr. In this case, the effective 
luminosity of an association of age t that will eventually form a total of A/"*h massive stars is 

M a ( \ 

£FW,H 2 (-A/"*fc,£) = \ /esc[* - (j ~ 1) tg] zl ^FUV,H 2 (^j) } • (17) 
j = 1 I *>*liveGM) J 

Here the first summation is over the M g < 5 generations or bursts that occurred during the 
association age t [i.e., those for which t — (j — l)t g > 0]. The second summation is over all the 
stars born in generation j and that are still on the main sequence at time t. (The stars in each 
generation are labeled in order of decreasing mass, and iu ve (j,t) labels the most massive star from 
generation j still on the main sequence at time t.) In practice, the second summation can be 
stopped at m < 2Mq due to the negligible FUV luminosities of such stars. 

Figure 4 shows the evolution of the FUV luminosity £fuv (N*h > t) for associations characterized 
by generations with A/"*/i/5 = 200, 20, and 2 high-mass stars (i.e., the associations ultimately 
form 1000, 100, and 10 high-mass stars, respectively). The decreases in £fuv are caused by 
the deaths of the most massive stars born in the early generations. The first generation is born 
at age t = 0. The initial increase is due to the gradual dispersion of the dust from the natal 
cloud around the first generation. The final decrease is due to the death of OB stars after the 
last generation has formed. The right vertical axis gives the contribution to the energy density 
£/fuv at a point located 100 pc away assuming that there is no intervening extinction [i.e., 
J7 FUV = (l/4vrr 2 c)£ F uv/(2070 A-912 A)]. Fi gure 4 shows that the maximum luminosity of an 
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association is attained at about the time the fifth generation dissipates the absorbing material 
around it (i.e., about 18 Myr after the first generation begins to form). At this time, the mean 
FUV luminosity per high-mass star formed (in the five generations) is L max f l ~ 1.6 x 1O 4 L0. 
Note that L maiX j L is only about 28% of the value (-Lfuv)/i calculated in §2.2 since many of the 
most luminous stars in the first four generations have died by this point. The contributions from 
generations 1 to 5 to L maXi/l are respectively 5.4%, 7.9%, 12.4% 21.6% and 52.7%. At 20 Myr, the 
mean FUV luminosity per high-mass star formed drops to ~ 1.0 x 10 4 L Q . 

3.2. Size Distribution of Associations 

Using compilations of radio HII regions in the Galaxy and optical studies of HII regions in 
nearby galaxies, MW97 showed that the luminosity distribution of giant radio HII regions in the 
Galaxy can be fit by a truncated power law of the form N a {> S) = N au [(S u / S) — 1], where S is 
the initial ionizing photon luminosity, N a {> S) is the number of associations with luminosity at 
least S, S u is the upper limit of the distribution, and N au is the number of the most luminous 
associations, with luminosities between 0.5 S u and S u . 

The luminosity distribution of OB associations does not continue as a power law in 
S to low luminosities because when the luminosity is dominated by a single star, the 
distribution is determined by the IMF. MW97 assumed that the underlying dependence on the 
number of high-mass stars in the association did not change at this point, however, so that 
M a {> N*h) = Nau[{N*h,u/ N*h) — 1]- Here N*h is the number of high-mass stars in all the 
generations of the association, and J\f*h,u is the number of high-mass stars in an association in 
which each generation has an ionizing luminosity S u (MW97 argued that the different generations 
of an association had comparable ionizing luminosities). If the lifetime of the association is t a , the 
birthrate of associations in which at least A/"*fc high-mass stars form is then 

-lf$5H- (18) 

MW97 assumed that associations had 5 generations, each lasting 3.7 Myr, so that t a = 18.5 
Myr; we shall adopt this value here. For the Galaxy, MW97 estimate an upper bound 
S u = 490 x 10 49 photons s _1 , which corresponds to S u /(s)h = 1180 massive stars per generation, 
or Af*h,u = 5900 massive stars over the life of the association. This is somewhat smaller than the 
value Af*h,u = 7200 estimated by MW97 because of the flatter IMF we have adopted here. The 
lower limit of the distribution is N*h,l = AC/i,u exp[— St/(No,u Su)], which is obtained by requiring 
that the integrated ionizing photon emission by the entire distribution of associations match the 
observed value. With S T = 2.6 x 10 53 photons s" 1 and N au = 6.1 (MW97), we find N* h ,l - 1-0. 

The birthrate of associations per unit disk area is <^ a (> AC/i) = A/" a (> N*h)/A c s. At the solar 
circle the effective area is A e g ~ 530 kpc 2 , so that the birthrate of associations per unit disk area is 

a> Mh) = 6.2 x 10- 4 (™ - l) kpc" 2 Myr" 1 . (19) 
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3.3. The Typical Star Formation Rate is Below Average 

Equation (|19|) reveals that, on average, one association (of any size, N*h,l < -^*h < ^~*h,u) 
forms in one kpc 2 each ~ 0.3Myr at the solar circle. The number of massive stars formed in 
a given area A and over a given time interval At will depend on the size N*h and number M a 
of associations formed. For small areas (e.g., ~ lkpc 2 ) large fluctuations are expected even if 
At is tens of millions of years. We shall show that as a result, the typical rate of massive star 
formation ?*/i,typ can be significantly less than the average rate in the Galaxy for small values of 
AAt/(A eS t a ). Since the FUV-ISRF is determined by the average rate of massive star formation in 
a small fraction of the disk over a relatively short time, this means that the typical value of the 
radiation field is significantly less than the average value. 

The number of associations M a (> N*h) with at least N*h high-mass stars that are expected 
to form during a period of time At inside an area A is 

AM) - *- (> ?- tMA ' ■ (20) 

Conversely, the largest M a associations formed in an area A and time At, can be expected to have 
at least . 

= - N A -"-; (21) 

high-mass stars, where we have used equation (^) to evaluate N a . For example, consider a large 
sample of sets of associations formed in AAt; the median size of the largest association in each set 
will be M*h(Na = ^) (i.e., it is expected that half of the sets of associations formed would have an 
association with N*h larger than this value). As a numerical example, A = 1 kpc 2 and At = 25 
Myr, give N*h{N a = 5) = 178, which is substantially less than N*h,u = 5900. 

This analysis carries over into comparing the mean and median rates of star formation in a 
given area of the Galaxy. The important point is that, averaged over a finite region of space and 
time, the star formation rate is most likely a relatively low value, but that if one waits long enough 
or chooses a large enough area, a very massive association will form in this region that raises the 
mean rate above the median or typical rate. For an area A and time interval At, the typical star 
formation rate can be estimated as the rate of formation in all but the largest association plus the 
rate of formation in the typical largest association, 

i*h - -7- K h dN a + -^-K h {\)N a [K h {i)]. (22) 

The ratio of this typical rate to the mean rate can be evaluated with the aid of equation (|T^), 
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S*/i,typ 



■Nmh.1 
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A/ *h u 



11 \ ACfoj 



(23) 
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For example, if the typical star formation rate is measured by observing the O and B stars, 
then we are only probing a time At ~ 25 Myr, and the typical to mean ratio within 1 kpc is ~ 0.86. 
We show in §4 that most of the typical FUV-ISRF originates from associations closer than ~ 400 
pc and younger than ~ 30 Myr; in this case, the ratio of typical to mean is 0.69. Moreover, when 
averaging over very long periods of time, about half of the local mean FUV-ISRF is produced by 
very infrequent large associations within less than H* ~ 85 pc, the observed scale height of OB 
stars. The ratio of the typical to mean star formation within H* during a period of 30 Myr is 
~ 0.34. In other words, for very long periods of time there are no large associations within H* and 
the FUV field is relatively low. Then, sporadically, a large association forms within and raises 
the field so much that, even though its duration is brief, it pushes the mean significantly above the 
median. Properly accounting for this effect, the typical FUV-ISRF is about half the mean value! 

As a sidenote to this section, it is interesting to compare the predictions of the MW97 model 
with observations of nearby OB associations. In total, about 7 associations younger than 25 Myr 
are observed within 1 kpc from the Sun. Equation ( |20| ) predicts that the 7 associations will each 
produce more than ~ 40 high-mass stars. The completeness of the observed sample is unknown, 
but the observations appear consistent with the model. MW97 came to the same conclusion 
by a slightly different approach: They calculated that there should be 7.3 associations within 
a kpc of the Sun that currently have at least 30 high-mass stars (i.e., excluding some young 
associations that have yet to reach their full complement of high-mass stars) and are younger 
than 25 Myr. In summary, the MW97 model for the OB associations is in good agreement with 
the observed associations in the Solar neighborhood. However, the observational bias is such that 
these associations are large and young. A significant contribution to the FUV-ISRF comes from 
smaller and/or older associations that are less readily detected. We shall use the MW97 model to 
describe these associations as well. 

4. Time dependent FUV-ISRF 

To calculate the evolution of the FUV-ISRF at a given point in the Galaxy, we assume that 
the birthrate of associations per unit volume has a Gaussian distribution with respect to the 
height z above the galactic plane. We take the scale height of the birthrate to be = 85 pc, 
the observed scale height of molecular gas and OB stars (Ferriere 1998, and references therein). 
Since the surface number density of young associations is small, the above distributions must 
be regarded as probabilistic. In the plane the distribution is assumed to be uniform (random) 
since almost all the FUV-ISRF originates within ~1 kpc, the scale length for dust attenuation of 
FUV, and therefore radial gradients in the birthrates are negligible. Large scale patterns of star 
formation, such as spiral arms, are present in the Galaxy, but it is not necessary to include the 
corresponding spatial variation in the star formation rate since the FUV field at a given point is 
mainly determined by the population of massive stars in a rather small region around it (within ~ 
0.4 kpc). Our adopted star formation rate is based on the global model of MW97 and is therefore 
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an average over arm and interarm regions of the Galaxy. Since the fluctuations in the ISRF we 
calculate are purely statistical, they are a lower limit to the fluctuations that occur as a point 
moves between arm and interarm regions. 

Figure 5 shows the projection on the plane of the positions of associations formed in a period 
of 500 Myr (left panel) and 30 Myr (right panel). The birthrate corresponds to the average value 
at the solar circle. Note that although there are more small associations than large associations, 
there are equal total numbers of stars in small and large associations. For a point at the center, 
most of the FUV radiation comes from within the circle. The FUV flux varies as function of 
position (or time), depending on whether the given point is close to large associations. 

The distribution of attenuating dust with respect to z is assumed to follow the distribution of 
HI as given by Dickey & Lockman (1990), i.e., 

n H (^) = 0.566 [0.69 e- (2/127pc)2 +0.189 e~ (z/318pc)2 +0.113 e- |2|/403pc ]cm~ 3 . (24) 

Note that bh is the space-averaged number of H atoms in a large volume that might contain a 
significant amount of hot gas (McKee & Ostriker 1977). The molecular gas is not considered 
due to its small filling factor. The dust optical depth to a given association is then assumed to 
be proportional to the column density of H atoms, iV(HI), along the line of sight. We adopt 
the physical properties of dust proposed by Witt and Gordon (1999) to estimate the average 
absorption and scattering by dust grains in the 912 - 1100 Aband and in the 1100 - 2070 Aband. 
The scattering and absorption optical depths are given by t S)CL = K Sja iV(HI), where the values of 
K S)a are given in Table Bl. The assumption of a constant extinction per unit column density is 
justified in statistical studies of the expected radiation field, but this approximation is not useful in 
calculating the contribution of a particular source because the extinction per unit column density 
is known to depend on the particular line of sight. Scattering removes photons from the direct 
line of sight to the associations, but adds photons from other directions. The effect of scattering is 
considered in Appendix B. 

The contribution of an association labeled 'V to the radiation density at a distance r can be 
approximated as 

n band « = r V%^e— W , (25) 
4irr 2 c A band 

where c is the light speed, and £band(0 is the effective luminosity (eq. |17|) of the association 
labeled i, which has ]\f*h high-mass stars and whose age is t. The term e _Tband accounts for both 
absorption and scattering (see Appendix B). 

At a given time t s i m of the simulation, the total energy density in each band ?7band(*sim) is the 
sum of the contributions MbandW from all the associations formed between t = and t = i s i m . Due 
to the 1/r 2 term in equation (|25|), the resulting radiation density [/band at a given point is very 
sensitive to the spatial distribution of the associations; in particular, a nearby luminous source 
can dominate the radiation field. The effect of extinction is very important for distant sources. 
Without extinction the contribution of FUV sources in a disk grows logarithmically with distance; 
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with extinction the total FUV field comes from sources within a few optical depths. In practice it 
is sufficient to consider only associations within about 1 kpc for conditions at the solar circle of 
the Galaxy. Additionally, old associations (with ages exceeding the lifetime of stars with masses 
~ 2 — 3 Mq) make a negligible contribution to J7fuv- These two facts are useful in reducing the 
number of associations that must be tracked in simulations to obtain an accurate accounting of 
the time dependence of the radiation field at a given point. 



4.1. The Mean FUV-ISRF 

Before presenting the results from simulations, it is useful to estimate the mean FUV energy 
density assuming a smooth emissivity per unit volume corresponding to the PDMF and a smooth 
spatial distribution; hereafter the "continuous method". This approach is necessary because the 
"discrete method" used in the simulation (assuming that FUV sources are highly discrete in 
space and time) is useful to get the median value of the radiation field, but is very inefficient 
in calculating the mean value because of the large integration times needed to obtain a good 
sample of rare events. However, the "continuous method" does not show the interesting time 
dependence of the FUV-ISRF. Knowledge of the ratio between the median and the mean value 
of the radiation field is important because low resolution observations (for example, in external 
galaxies) that average over huge spatial volumes can provide an estimate of the mean FUV-ISRF, 
but the interstellar gas is subject typically to the median FUV-ISRF. 

The mean energy density < ?7band > in the plane can be calculated by adding the diluted 
extincted and scattered radiation from a smooth distribution of stars following a PDMF. Let 

d 2 N*(m,r) 

n *( m > r ) = -H7T, ( 26 ) 

aVdlnm 

be the PDMF per unit volume at a point r. The mean energy density of radiation at the origin is 
then 

\ f g — Tband(r) rtUu 

(t^band) = — r / — j — o~dV / n*(m,r)L band (m)(ilnm , (27) 

cA band J 4irr z J„ H 

where Lband (jn) is the mean luminosity of a star of mass m after subtracting the fraction lost in 
the parent cloud. As shown in Appendix B, we approximate the optical depth, including both 
absorption and scattering, as Tband = T a + r s (l — 1Z), where 1Z is a factor that is somewhat less 
than unity for scattering that is strongly peaked in the forward direction. 

Since the stars that dominate the FUV flux are relatively nearby, we can ignore the Galactic 
radial gradient in the density of stars and express the volume density of stars in terms of the 
surface density as 

ra*(m, r) = s*(m) D(z) (28) 
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where z is the distance normal to the plane and, for a Gaussian distribution, 

D(z) 



(29) 



Recall that for a constant birthrate, the PDMF for stars younger than the age of the disk is related 
to the star formation rate by s*(m) = ?*(m)t ms (m). Equation (p7|) can then be rewritten as 



(0 



4VT Jband S Li- 



band/ 



cA 



band 



where 



c(m) L haxiA (m) dhim 



(30) 



(31) 



is the mean rate of emission of radiation (after subtracting the absorbed emission in the natal 
clouds) per unit disk area in the band, and 
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(32) 



is 4-7T times the mean intensity due to unit FUV emissivity in the disk. 

At the solar circle S^ H2 ~ 2.5L0/pc 2 and £l fuv ~ lO.9L /pc 2 . Numerical integration 
of equation (32), taking into account multiple scattering (Appendix B), gives respectively 



0.88 and 1.19 for the value of 47r l 7band in the H2 and FUV— H2 bands. From equation (30) 
one obtains {U^z) = 15.7 x 10~ 17 erg cm~ 3 A \ (£/fuv-H2) = 13.8 x 10~ 17 ergcm -3 A , and 



(C/fuv) 
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Hereafter, we denote these mean values as (?7band) 
density for an infinite period of time. 
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since they correspond to the mean energy 



To clarify how different regions contribute to the mean radiation field, consider a simple model 
consisting of two regions: the first is a sphere of radius centered at a point in the midplane, 
and the second is a zero thickness disk that covers the midplane but excludes a radius around 
the point in question. The volume number density per unit mass of FUV sources in the sphere is 
n*(m) = s*(m)/2H*, and the surface density in the disk is s*(m). The UV extinction is assumed 
to be e~ kr , where k = [n a + k s (1 — TV)]n^ ~ 7 x 10 -22 cm -1 (see Appendix B). Then, 
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(33) 



rn, 



The two terms in the bracket correspond to Air J in this approximation. Since kH* <C 1, the 
contribution of the sources in the sphere to the function Air J is ~ 1/2. This term represents 
sources out to a distance of order the scale height H* of the disk. Beyond H*, the contribution of 
the sources in the disk is \ Ei(kH*), where E\ is the first exponential integral. Since kH* ~ 0.18 
and -Ei (0.18) ~ 1.3, the second term is about 0.65. In other words, in calculating the mean field, 
we find that roughly half the radiation comes from sources closer than ~ = 85 pc, and half the 
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radiation arises from more distant sources in the disk. However, as shown in equation (23), the 
typical number of FUV sources inside of 85 pc is much less than the mean value. Therefore, most 
of the time, the FUV radiation field is expected to be well below its mean value. We also see that 
beyond r > k^ 1 ~ 0.46 kpc the contribution to the FUV field rapidly falls off due to attenuation 
of the distant FUV flux by dust. 

In order to quantify the simulation time necessary to get a good estimate of (U), we calculate 
the accumulated typical FUV radiation field during a period of time t smi , 
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where the typical maximum value of A/*^ in the area A during time interval At = t s im is given by 
equation (|2l|). This equation applies only if At s i m is large enough that Af*h, typ max > ^V*/i,z> for 
smaller values of ^4i s im, f*typ(^) A, t S i m ) = 0. 

Note that in computing (J7band)t sim fo r a given spatial point, the area of simulation A s i m must 
be chosen so that Tb an d S> 1 at the outer boundary r max (see eq. (|3^)). Otherwise, the simulation 
area would not include significant contributors to the UV flux at the central point. This is 
equivalent to r max 3> k~ x , which for our standard values of «h and n a ^ s results in r max 3> 0.46 kpc. 
Assuming a square simulation area, this corresponds to A s \ ra S> 0.85 kpc 2 . Given that A s \ m ~ 5 
kpc 2 , the simulation area is no longer a factor in (U); only the simulation time enters. Figure 
6 shows the ratio ( CAband )t sim /( CAband )oo as function of the time i s i m . Extremely long times are 
required to get a good estimate of (?7band) from the simulations. 



4.2. Simulation in the Galactic Plane at the Solar Circle 

Over a sufficiently large area ^4 s i m to include all possible contributions to the FUV field not 
severely attenuated by dust, the associations form following the truncated birthrate law given 
by equation (|^). We assume a constant interval of time At a form between the formation of 
associations in the simulation area A^im- Equation ( |i"9|) gives this interval as 

At aj form = ~T~N \ ' 

Mau (^Sf -1 ) Asim 
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Substituting M au = 6.1, t a = 18.5 Myr, ftf*h,u = 5900 and N*h,l = 1 ( see §3.2), we obtain 
At a = 514 (A e g/A s i m ) yrs. At the solar circle A e g = 530 kpc 2 and we simulate a disk area 
^sim = 5 kpc 2 , so that an association forms in A s - im every 5.5 x 10 4 years. 

In the simulation, the position of each association is taken to be completely random in the 
plane (x,y), but in the vertical direction (z) the random distribution is weighted by a Gaussian 
profile of scale height H* = 85 pc. The masses of the stars in each generation of each association 
are chosen using equations (Al) and (A2). The adopted gas distribution from which we derive the 
dust distribution is given by equation (|2^), and the methods for including dust absorption and 
scattering are given in Appendix B. To avoid boundary effects, the evolution of [/band is studied 
only in the central part of ^4 s im- 

Figure 7 shows the time evolution of C/fw a ^ a P orn t located in the center of simulation 
area and in the galactic plane. The radiation field energy density U-n is expressed in units 
of 10 _17 erg cnr 3 A . We note that Habing (1968) found that the local interstellar radiation 
field is approximately t/jPuv = ^-6 X 10~ 17 ergcm~ 3 A . The first association in the simulation 
volume forms at i s i m = 0. Although difficult to see in Figure 7, it takes ~ 40 Myr to build up 
the FUV-ISRF once star formation commences. After this time, a "steady state" is achieved 
with large time fluctuations of U-n- In addition, large spatial gradients are derived when U-n 
is evaluated at different points in the plane at a given instant of time. Note that although the 
probability is small, a massive star can form in a small association (see Appendix A). These small 
associations are abundant, and many of the HII regions encompassing the point are produced by 
these unusually massive stars in small associations. 

The differences in the behavior of the radiation field in the H2 band [912 — 1100 A] and in the 
FUV— H2 band [1100 — 2070 A] are shown in Figures 8a,b. Note the high concentration of points 
around the median value, which corresponds to Uu2 ~ f T FUV-H2- The ordered sequence of points 
outside the crowded region in Figure 8b corresponds to periods of time when the radiation field is 
completely dominated by one association. In Figure 8b, two larger excursions to high values of U 
have been connected by lines in order to show how the field evolves as generations form and stars 
die. The large fluctuation that starts at t s i m = 1875 Myr is produced by an association at 190 
pc with 2830 high-mass stars. The fluctuation that starts at t s { m = 3975 Myr is produced by an 
association at 32 pc with 440 high-mass stars. The main trend seen in Figure 8b is that, as C/fuv 
increases, the ratio Ua 2 to C/fuv-h 2 increases. This expected correlation arises because high Z7fuv 
is often associated with massive (and therefore hot) nearby stars. 

Figure 8c shows the fraction of £4,and in the [1100 — 2070A] band that comes to the observer as 
scattered light (see Appendix B), as function of U-n- The median value of f/scatt/^tot is ~ 0.18. 
We discuss this median value in comparison with observational data in §5. 

Figure 8d shows for the FUV band the fraction of total field U-n that comes to the observer 
from the source that provides the larger contribution. For values of U-n around the median (~ 7) 
the dominant source provides in average 20% of the energy density U-n, but the values range 
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from 10% to 50%. For large values of U—n the radiation field is totally dominated by only one 
source. The distribution of points in Figure 8d indicate the degree of anisotropy of the field as 
function of U-n- 

In order to calculate the mean and median values of the FUV radiation field in the plane, a 
very long simulation was performed with t s j m ~ 57 Gyr. Excluding the first 400 Myrs, Figure 9a 
shows the corresponding distribution of the "fraction of occurrence" of [/band va l ues - The fraction 
of occurence is defined as —dP(> U)/dlnU, where P(> U) is the fraction of the simulation 
time that the point has an FUV energy density exceeding U. The fraction of occurrence can be 
thought of approximately as the probability of lying within a factor of 1.6 of [/band- The arrows 
indicate the mean and median values of the distributions. The two arrows for each mean value 
indicate the value including (the higher) or omitting (the lower) the time steps when the point is 
inside an HII region. The difference with and without HII regions is very small for the median 
because a given point is in an HII region only about 2% of the time. The difference for the 
mean is larger because [/band * s vei T high when the point is in a HII region. The expected mean 



value of [/band (including HII regions) from equation 34 and Figure 6 at 57 Gyr of simulation 
is ~ 0.9 of the theoretical asymptotic value calculated in §4.1. Therefore, even when 57 Gyr 
of simulation are considered, the simulation does not provide an accurate estimate of the mean 
values because extremely rare events with large associations very close to the point in question 
tend to be missing. On the other hand, the mean value of [/band excluding HII regions is closer 
to its asymptotic value since regions of space very near the association are excluded by the 
requirement that the point be outside the HII region. Note that the energy density in the FUV 
band is never lower than {Jftjv mm = 2.4 x 10~ 17 ergcm~ 3 A ; in the H2 band the minimum 
value is 1.5 x 10~ 17 erg cm -3 A . These minimum values correspond to the periods of time when 
no young associations are present in the central region of the simulation volume. 

Figure 9b shows the probability P(> U) that, at a given instant in the simulation, the FUV 
energy density is greater than U. As stated above, the energy density in the FUV band is never 
lower than [/fuv, min in the simulation, and therefore P(> C/fuv, mm) = 1- The curves labeled 
1 Gyr, 4 Gyr, and 56 Gyr, correspond to the probabilities P(> U) calculated from simulations 
extending over time intervals At = 1 Gyr, 4 Gyr, and 56 Gyr, respectively. Note that the 
peak energy density decreases as At drops. Even for a simulation time of 56 Gyr, events with 
U-n > 4000 are missing. For U-n > 10, the probability P(> U) can be fitted by the power 
law P{> U) ~ IQUZn 2 . For E7_ 17 > 300 a better fit is P{> U) ~ 15 EC?/ 3 . This dependence 
arises from the fact that the volume of the sphere around a source of luminosity C where the 
energy density exceeds U is proportional to (£/[/) 3 / 2 . The volume filling factor </>(> U) of these 
spheres therefore scales as [Z~ 3 / 2 . For <p « 1, we have P(> U) = 4>{> U), so that in this case 
P{> U) oc [Z~ 3 / 2 . This relation is not exact because, for example, extinction is neglected, and the 
overlapping of spheres from different sources is not considered. However, the fit becomes better 
at high values of U since these effects become unimportant there, since the volume around the 
sources is then small and consequently overlapping and extinction can be neglected. Figure 9b also 
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shows the typical time t\j (right vertical axis) between two events that raise the energy density to 
> U at a fixed point in space. It can be shown (Hollenbach et al. 2002) that tjj = tp\jy/(j)(> U), 
where iFUV is the typical time that a source (an association) emits strongly in the FUV band (see 

3 11 

Fig. 4). The power law tjj = iFixv/-P(> U) = 2U_ 17 Myr assumes £fuv = 20 Myr and is good for 
U- 17 > 10. For U- 17 > 300 a better fit is t v = 1.3 UlQ Myr. From this relation one can see that 
for t sim = 57 Gyr only the events with U-n < 1000 are well sampled. 

For the FUV band, Figure 10 shows the evolution of the mean (U) and the median C/ me d 
radiation fields (including HII regions) normalized to the theoretical value (U)^ derived in the 
previous section. The calculation of these averages is over the period of time t s i m starting at time 
zero when associations begin to form in the simulation area. The theoretical ratio {U)t sim / {U)^ 
(plotted in Fig. 6) is also shown in Figure 10 for comparison. Note that the median initially grows 
for ~ 200 Myr and then stabilizes, showing that the simulation is very efficient for calculating the 
median value of the field. The amplitude of fluctuations tends to decrease with time since the 
averages are made over the entire history of the simulation. The large fluctuation in (C/sim)/(C^)oo 
at ~ 100 Myr is due to an association of ~ 300 high-mass stars (the most massive star with a 
mass of 86 M ) born at t s j m ~ 60 Myr at a distance of ~ 120 pc (with z ~ —80 pc) from the 
central point. At t s i m ~ 75 Myr, U ~ 25 x 10~ 17 erg cm~ 3 A \ mainly due to the FUV flux from 
this association. The second large fluctuation at ~ 600 Myr is due to an association of ~ 41 
high-mass stars born at 

^sim ~ 540 Myr at a distance of ~ 30 pc (with z ~ +23 pc) from the 
central point. This association is relatively small, but its most massive star has a mass of 114M Q . 
At t sim ~ 560 Myr, U ~ 80 x 10~ 17 ergcnr 3 A" 1 , mainly due to the FUV flux from this star. The 
large fluctuation at ~ 21 Gyr is due to a very rare event: an association of 3120 high-mass stars 
at a distance of ~ 35 pc (with z ~ —9 pc) that raises the FUV radiation field to a peak value of 
U-rj ~ 3200. The radiation field has U-n > 1000 for over 20 Myr, enough to raise (f7 s i m ) over 
{U}oo even when the average is made over 21 Gyr. 

As stated in §3.3, the star formation rate is expected to fluctuate during the simulation. 
Figure 11 shows the evolution of q^, the surface (number) density of high-mass stars in the 
simulation. The thick line corresponds to the average surface density over A sim (i.e. 5 kpc 2 ), 
whereas the thin line corresponds to the average within 500 pc from the central point. As expected, 
the fluctuations are larger for the surface density within 500 pc. For a simulation interval between 
0.1 to 57 Gyr, the mean value of over j4 s im is ~ 540 kpc 2 , and within 500 pc is ~ 530 kpc 2 . 
Due to the large simulation time, these mean values of are similar and correspond to a SN II 
every 62 yrs in the whole galaxy (we take A e R = 530 kpc 2 ). 



5. The Local FUV Radiation Field: Observations and Theory 
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5.1. The Diffuse FUV Radiation Field 

Observationally, it has been established that diffuse radiation coming from regions where no 
sources are observed (i.e. far from the regions of recent star formation) is a minor component 
of the FUV-ISRF. This diffuse component arises primarily from large angle scattering by dust 
grains of FUV photons from Galactic sources, with a smaller contribution from extragalactic 
sources. Measurements of the diffuse FUV radiation field from these regions place an upper 
limit of 200 - 400 photons s" 1 sr" 1 in the [912 - 1200 A] band (Murthy et al. 1991; 
Bowyer 1991; Henry et al. 1980). This upper limit corresponds to a radiation density 
^diffuse ~ 1-7 — 3.4 x 10~ 18 ergcm~ 3 A and represents less than 10% of the estimates of local 
FUV energy density. By comparing a map resulting from summing the fluxes from individual 
FUV bright stars (Gondhalekar, et al. 1980) to a map constructed from large angle detections, 
Henry (1991) demonstrated that the FUV radiation field is dominated by the direct (or forward 
scattered) flux from relatively nearby FUV-producing stars. In the previous section we presented 
model results that showed that the median value of ?7 sca tt/C^tot is ~ 0.18, which is higher than 
the above estimates of [/diffuse /t^tot ^ 0-1- This difference arises because the "diffuse field" 
[^diffuse comes from regions in the sky far from strong stellar sources, whereas the scattered field 
U scs ,tt refers to all scattered flux, including that scattered through small angles. In fact, we find 
(Appendix B) that about 1/3 of the once-scattered radiation from a source comes from lines of 
sight within a cone of five degrees from the observer to the source. Taking this into account, we 
conclude that our estimate of the scattered flux is consistent with the observed diffuse field. 



5.2. The Observed FUV Radiation Field in the Solar Vicinity 

Neglecting diffuse radiation, the spectral energy density U\ at wavelength A in the immediate 

solar vicinity can be directly measured by adding the observed flux F\(i) of all the FUV stellar 

sources i. By including all the data in the S2/68 Survey (more than 50,000 stars), Gondhalekar et 

al. (1980) estimated the flux at several wavelengths in the FUV band. At 1400A the estimated 

energy density is about 5 x 10 _17 erg cm~ 3 A 1 . Habing (1968), like Gondhalekar, added the fluxes 

from individual stars to estimate U\ in the immediate solar neighborhood. The estimated value for 

U ° was 4.2 x 10 _17 ergcm~ 3 A when the direct light from all the bright stars B9 and earlier 
1400 A to to 

were included. For the extreme case of albedo 0.9, the Habing's value with diffuse light included 
was C r 14Q0 ^ = 7.5 x 10 _17 ergcm~ 3 A . However, this value depended on the assumptions made 
to infer the FUV flux for each object from its visual magnitude my and its color excess E(B — V). 
Jura (1974) recomputed U\ near 1000 A and obtained a value ~ 7 x 10™ 17 erg cm~ 3 A 1 by using 
updated calculations of stellar atmospheres and including all stars B5 and earlier listed in the 
Catalogue of Bright Stars. Utilizing a different technique, Henry et al. (1980) have directly 
measured the local U\ by taking large beam ( 12° x 12°) observations over about one third of 
the sky and extrapolating to obtain the unobserved portion. Henry et al. found that the average 
energy density is 5.8 x 10~ 17 ergcm~ 3 A 1 in the band [1180—1680 A], and 5.3 x 10~ 17 ergcm~ 3 A 1 
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at 1400 A. Draine (1978) compiled a number of the observations and models of that epoch, and 
did a fit to some of them (primarily Henry et al (1977), Witt and Johnson (1973), Belyaev et 
al (1971), and Jura (1974)) to obtain what has become another standard (besides Habing) field, 
sometimes called the "Draine field". Averaged over the FUV band, it is approximately 1.7 times 
higher than the Habing field, i.e., the Draine field has Go = 1.7. A summary of these estimates is 
given in Table 2. 

The asymmetry of nearby FUV point sources results in a non- isotropic local FUV field. Most, 
if not all, OB stars are born in associations, and there are few associations that are young and 
close enough to produce a significant contribution to the FUV field. A sample of the nearby OB 
associations reveal about 7 associations within a kiloparsec of the Sun (see §3.3). The distribution 
of their distances above and below the galactic plane is quite symmetric, but their distribution in 
the plane is totally asymmetric. When the observed FUV field is plotted in Gould coordinates 
(Galactic Coordinates tipped by about 19°), a clear excess toward the active section of the Gould 
Belt is apparent (Henry et al. 1980). Even when our simulations do not consider correlation 
between the positions where coeval associations born, the results in Fig 8d shows that because 
FUV sources form in groups, FUV field is in general non-isotropic, in particular when the radiation 
field is high. 

5.3. The Mean and Median FUV Radiation Field at the Solar Circle 

Complementing these more or less direct determinations of U\ in the immediate solar 
neighborhood are the estimates, such as ours, obtained by modeling a spatial distribution of the 
relevant spectral type stars and the dust. The advantages of modeling are that it corrects for the 
incompleteness in the catalogs, it is useful to quantify the variations from different distributions, 
and it generalizes to other locations than the immediate solar neighborhood. The disadvantage 
is the necessity to introduce simplifying assumptions such as smooth dust distributions and 
prototype stellar luminosities. However, the "direct" procedure is only useful to give a snapshot 
of the local FUV-ISRF at the current time, whereas the theoretical models permit the study of 
the temporal and spatial behavior of the UV field in the interstellar medium, which is the main 
subject of the present paper. 

In the previous section, the mean FUV field was calculated assuming the model ingredients, 
but, for the sake of comparison, we reproduce here one of the classical calculations of Habing. 
The "continuous approximation" in the Habing (1968) model includes the contribution of 10 
associations and 577 field stars between BO and B5. The FUV emission of these sources is assumed 
to be uniformly distributed in the plane (the field stars in 1 kpc 2 and the associations within 
1 kpc). The distribution is proportional to exp[— (z/120 pc) 2 ] in the vertical direction. Habing 
assumed a homogeneous plane-parallel layer of dust of thickness 150 pc producing an extinction 
of 3.5 mag kpc -1 . By integrating over this continuous medium (there are neither stars nor 
associations, but an equivalent emitting medium), he obtained what we term the observed-source 
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mean radiation energy density at 1400A of (^ 1400 ^) = 3.4 x 10~ 17 ergcm~ 3 A -1 . (This value 
assumes that the dust has zero albedo, so there is no scattered light. For an albedo of 0.9, Habing 
obtain (U °) = 5.6 x 10 _17 erg cm -3 A" 1 ). The distinction between Habing's observed-source 
mean and our mean is that ours is a mean based on the star formation rate in the entire Galaxy, 
whereas Habing's is based on the observed FUV sources in the solar vicinity. The surface 
emissivity £l fuv in Habing's model is about 1/3 of the value in our model; this is plausible, since 
we expect his value to be closer to our typical value than to our mean (§3.3). Because his dust 
opacity is somewhat higher than ours, his observed-source mean radiation density is smaller yet, 
about 1/4 of the value of (Cfuv) calculated in §4.1. 

The observed-source mean FUV energy density estimated by Habing (1968) and the mean 
FUV energy density calculated in this paper are summarized in Table 3, while the typical values 
are summarized in Table 4. Habing's typical value is different from ours: his estimate is obtained 
by excluding the space in the neighborhood of associations (about 10% of the total volume) and 
taking the average over the remaining space. As discussed in §4, we estimate the typical FUV 
radiation by calculating median value of the FUV field at a point in a large simulation of the 
model. Compared to the band-averaged FUV "Habing field", our model results give Go = (1.6,3.1) 
for the (median, mean) radiation field at the solar circle. Table 4 also gives the value below 
and above which the FUV energy density exists 20% of the simulation time. A more complete 
characterization of the FUV radiation field is given by the probability distribution in Figure 9. 
This representation obscures the temporal correlations of the radiation field, however, so it is also 
useful to use the FUV field history shown in Figures 7 and 8a to characterize the typical behavior 
of the FUV field. 

6. Effect of Runaway OB Stars and Supernova 

6.1. Runaway OB Stars 

"Runaway" OB stars, as defined by Blaauw (1961), are OB stars with peculiar space velocities 
that exceed 30 km s _1 and that are directed away from known clusters and associations. About 
10-30% of O stars and 5-10% of B stars are observed to be runaways. These runaway stars have 
large space velocities (30-200 km s _1 ), and most of them appear to be single (Gies 1987; Stone 
1991; Hoogerwerf et al. 2000a, b). The most viable scenarios for the formation of runaway stars 
are: (i) a supernova explosion in a massive binary (Zwicky 1957; Blaauw 1961), in which the 
primary becomes a compact object and the secondary moves away with a velocity comparable to 
its pre-explosion orbital velocity; and (ii) dynamical ejection from a dense stellar cluster (Poveda 
et al. 1967; Gies Sz Bolton 1986), in which the most efficient interaction is the encounter of 
two hard binary systems. Recently, milli-arcsecond astrometry provided by Hipparcos and radio 
observations have provided specific evidence for both scenarios (Hoogerwerf et al. 2000a,b), but 
the relative importance of the two mechanisms remains unknown. 



-26- 



Since the OB stars produce the FUV radiation field, and a non-negligible fraction of them 
are runaways, we must estimate their effect on the evolution of the FUV-ISRF. We assume that 
about 20% of O stars and 7% of B stars in each generation are ejected in random directions with 
a velocity weighted by the probability distribution 



with v c = a = 30 km s , appropriate for runaways. We consider only stars with masses larger 
than 8 M Q ; i.e. O and early B stars. Our primary assumption is that OB-star runaways (regardless 
of the ejection mechanism) are ejected at the moment they form. In the dynamical ejection 
scenario, the runaways are ejected very early in the cluster evolution, so that this assumption is 
well justified. In the binary-supernova scenario, the creation of the runaway is delayed until the 
explosion of the more massive component of the binary system. Therefore, the instantaneous 
ejection assumption overestimates the size of the sphere of influence of the SN-ejected runaways 
and provides an upper limit on the effect of runaways. 

To estimate the effect of runaway stars on the FUV radiation field, we compute the difference 
A[/_i7 jmn of the radiation energy density produced by the stars with masses over 8 M in 
two cases, with and without runaways. Figure 12 shows the effects of including runaways in 
the simulation performed in §4.2; that is, the simulation parameters are exactly those used to 
generate the time evolution of the FUV-ISRF in Figures 7 and 8, so that in both simulations the 
associations are formed with the same size, and at the same place and time. Additionally, we have 
assumed periodic boundary conditions; runaway stars that leave the simulation area from one 
boundary are reinserted at the opposite boundary. 

The heavy line in Figure 12 shows AC/_i7 jrun (t s i m ) for the same time interval in Figure 8a. 
The dotted line shows the total energy density ?7fuv from the simulation in §4.2 (i.e. considering 
stars of all masses and assuming that they all remain at their birth places). On average, at a given 
moment of time, there are about 200 runaway stars in the simulation area (~ 40 O runaways and 
~ 160 early B runaways in 5 kpc 2 ), but very few of them make a passage close to the observer and 
create a spike in the radiation field. Most of these spikes are produced by nearby associations; 
in the plot, only one significant spike is produced by a runaway coming from a distant (844 pc) 
association. For FUV spikes characterized by peak strengths U-n, the duration of the spike is 
of order At v ~ 10 6 L^v^Ul _l( 2 years, where L5 — Lpjjy /TO 5 L is the stellar luminosity and 
vj = v/100 km s _1 is the runaway speed. Spikes that dominate the FUV field occur roughly once 
every ~ 100 Myr. The negative values of AU-n ;Tun are also produced by nearby associations, 
since on average their runaway stars are at greater distances from the point in question than is the 
natal association. 

The median value of AU-n jTUn in the simulation is 0.016 (this slightly positive value is due 
to the fact that it is assumed that the parent cloud does not absorb radiation from the runaways). 
The ratio of the fluence in the spikes to that in the rest of the radiation field is ~ 0.025. Therefore, 




(38) 
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the conclusions drawn in previous sections (assuming that all the stars remain at their place of 
birth) do not change significantly if runaways are considered. However, very sporadically, the 
FUV-ISRF is dominated by the radiation from a runaway. These dominant spikes tend to occur 
when C/fuv is high (i.e., there are nearby associations). 

6.2. Supernovae 

Supernovae produce very large FUV luminosities for very short periods of time. Core-collapse 
supernovae produce UV radiation by two separate phenomena. One is the FUV radiation burst 
due to shock breakout, and the second is the FUV emission associated with the visible supernova. 
Scalo et al. (2001) estimate that each process typically produces about 10 47 erg in UV photons, 
over time intervals of less than an hour and 2 — 3 months, respectively. The luminosities in these 
phases are very large compared with the mean FUV luminosity per main sequence high-mass star 
(i.e. (Lpuv)fe = 5.8 x 1O 4 L , see §2.2) However, the luminosity-weighted lifetime of main sequence 
high-mass stars (i.e. (£fuv) = 7.8 x 10 6 yr) is much longer than the short periods where SN are 
actively emitting FUV photons, and as a result the total FUV energy emitted per main sequence 
high-mass star ((Lfuv)/i x (*fuv) ~ 6 x 10 52 erg) is about five orders of magnitude larger than the 
FUV fluence per SN (~ few 10 47 erg); indeed, the emitted FUV energy is about 60 times larger 
than the typical total energy released by a SN (~ 10 51 erg). Therefore, SN FUV spikes make a 
negligible contribution to the median and mean values of the radiation field calculated in previous 
sections. 

However, as pointed out by Scalo et al. (2001), these FUV spikes may be relevant for 
atmospheric chemistry and for inducing mutations in living creatures. Scalo et al. (2001) calculate 
the frequency of events exceeding the solar flux in the 2000 — 3000A at 1 AU and conclude that 
such events have occurred several times during the history of the Earth. Scalo et al. (2001) 
have considered an homogeneous distribution of SN. If as in the present work, massive stars 
are assumed to form in associations, high U events tend to occur in groups because high mass 
members of a given association explode within a relatively short period of time and at similar 
distances. Then, the time intervals between events (in different groups) exceeding a given flux 
increase by a factor s*h/s a (> 1) = lnN*h,u = 8.7. In addition, we note that the spectrum of the 
shock breakout is dominated by extreme ultraviolet and soft X-ray emission (Matzner Sz McKee 
1999), so only a small fraction of the emission is in the FUV band. When this energetic radiation 
strikes the Earth's atmosphere, much of it will be reprocessed to lower energies. The analysis of 
this process and of the effects on atmospheric chemistry is beyond the scope of the present work. 
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7. Conclusions 

The presence of a time and space-dependent FUV field may have a major impact on the 
structure of the neutral interstellar medium. In particular, the main process of H2 and CO 
destruction is photodissociation by radiation in the H2 band (912-1100 A), and the primary 
heating process of the neutral interstellar gas is grain photoelectric heating by radiation in the 
FUV band (912-2070 A). The FUV controls the amount of interstellar material in the cold (CNM) 
and warm (WNM) neutral phases. 

In this paper, we have presented a study of the temporal behavior of the radiation field in 
these bands, as well as at 1400 A. To simulate the evolution of the radiation field at the solar 
circle of the Milky Way, we have adopted: a) a power-law initial mass function with an upper 
mass cutoff for high-mass stars in associations (§2.1); b) a stellar-mass dependence of the FUV 
and H2 band luminosities, and the main sequence lifetime for stars of various masses (§2.2); c) 
a stellar-mass dependence of the ionizing luminosity from high-mass stars (§2.3), which, coupled 
with radio observations of the total ionizing luminosity in the Galaxy, determines the rate of 
high-mass star formation (§2.4); d) a history of the star formation rate in an association (§3.1); 
e) a dependence of the birth-rate of associations on the number of high-mass stars J\f*h i n the 
association (§3.2); and f) a vertical distribution of associations and dust (§4). Our calculations 
are based on the distribution of OB associations in the Galaxy as determined by MW97 and are 
completely independent of observations of the local interstellar radiation field. 

The most important conclusions are the following: 

1) We have computed the typical (or median) value of the FUV energy density at the solar 
circle in various bands and at 1400 A. Our computed value of the median FUV energy density 
at 1400 A at the solar circle (~ 7.5 x 10~ 17 ergcm~ 3 A^ 1 ) compares reasonably well with the 
currently measured local value at 1400 A (5.3 x 10 _17 ergcm -3 A -1 ; Henry et al. 1980) and with 
Habing's and Draine's recommended values for the typical local interstellar medium (~ 5 x 10~ 17 
and ~ 9.8 x 10 _17 ergcm -3 A -1 , respectively). The observationally estimated local value at 1400 
A is near the most probable FUV energy density values in the simulation (Fig. 9). Therefore, the 
local neighborhood has neither an exceptional high or low number of OB associations. We find 
band-averaged median energy densities of (7.4, 7.1, 7.2) xl0~ 17 erg cm" " 3 A" 1 for the ([912-1100 
A], [1100-2070 A], [912-2070 A]) bands, respectively. Our value for the FUV band is 1.6 times the 
value recommended by Habing (Go = 1.6), which is close to Draine's (1978) suggested value of 
Go = 1.7. 

2) MW97 proposed a distribution for the birth rate of associations in terms of the number of 
high-mass stars in the association, dSf a /d\nftf*h oc l/JV"^, based on observations of giant radio 
HII regions. Here we have shown that this distribution produces a median value of the solar circle 
FUV-ISRF in agreement with the currently measured local value. This is consistent with the 
extension of the MW97 distribution to small associations, as they proposed. 
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3) We have included scattering and absorption by interstellar dust. The absorption by dust 
limits the region that contributes to the FUV field at a given point to distances less than about 
500 pc. The scattering produces a "diffuse" field that is observed to be less than 10% of the total 
FUV field. Our models are consistent with these observations. Thus, most of the FUV energy 
density at a given point at the solar circle comes from the sum of discrete sources. 

4) We have computed the average (or mean) of the FUV energy density in various bands 
and at 1400 A. The mean value is nearly twice the median value, due to infrequent spikes in the 
radiation field. The mean energy density at 1400 A is 14.4 erg cm" ~ 3 A- 1 . The band -averaged 
mean energy densities are U- 17 = (15.7, 13.8, 14.2) for the ([912-1100 A], [1100-2070 A], [912-2070 
A]) bands, respectively. 

5) The fluctuations in the H2 band [912-1100 A] are larger than those in the FUV band, and 
the ratio of the mean to median energy density in this band is larger than that in the FUV band. 
This is both because the extinction in the H2 band is larger than in the FUV band and because 
the sources contributing to the H2 band are rarer than the those contributing to the FUV band. 

6) The FUV field undergoes substantial fluctuations, with the larger ones being less frequent. 
Figures 7 and 8 show that at the solar circle the energy density of the radiation U fluctuates 
above the median value by a factor 2 to 3 every 100 — 200 Myr. For U >> £/ me d, the time interval 
tu between episodes with fields greater than U is roughly given by tu ~ 2 U% Myr or 2 ulQ 
Gyr (Figure 9). These infrequent epochs of high fields last of order 30 Myr, the lifetime of the 
nearby OB association that creates the high field. As we shall show in a subsequent paper, these 
fluctuations affect the relative fractions of cold neutral medium and warm neutral medium, and 
thereby create a fluctuating state of the ISM as a function of space and time. 

7) A point at the solar circle is inside the H II region of a nearby OB association about 2% of 
the time, based on the assumption that the size of the H II region is that of a Stromgren sphere in 
a medium with a density of n = 1 cm" 3 . 

8) In general, a significant fraction of the FUV radiation field comes from a single source. 
For values of U-n around the median (~ 7) the dominant source provides on average 20% of the 
energy density, but the values range from 10% to 50%. For large values of U-n, the radiation 
field is totally dominated by only one source. Therefore, the FUV field is asymmetric at a 
given point, the asymmetry growing for higher fields. Very low values of the radiation field are 
very infrequent — indeed, values below 1/3 of the median FUV field are totally absent from our 
simulations. 

9) OB runaways produce short-lived spikes in the radiation field. At a given moment of time 
at the solar circle, we estimate that there are about 40 runaway stars per kpc 2 (~ 8 O runaways 
and ~ 30 early B runaways), but very few of them make a passage close to the observer and 
create a spike in the radiation field. However, very sporadically (on the order of once every ~ 100 
Myr, see Figure 12), the FUV-ISRF is briefly dominated by the radiation from a runaway. These 
dominant spikes tend to occur when C/fuv is high (i.e., there are nearby associations). However, 
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the duration of the spikes is short (about one Myr) and the total fluence in the spikes is negligible 
in comparison with the total fluence from the stars in the associations. 

10) Supernovae raise the radiation field significantly above the median value every ~ 5 x 10 4 
years. However, their duration is so brief that the total fluence in the spikes is negligible in 
comparison with the total fluence from the stars in the associations. 
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TABLE 1 
Stellar Properties 



Quantity Units Approximation Mass Range 
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TABLE 2 



The Local Energy Density 



Authors 


Method 


Wavelength 


U [10~ 17 erg cm -3 A -1 ] 


Hahine- 1 9fi8 


A 


1000A 


3.3 






1400A 


4.2 






2200A 


2.2 


Jura 1974 


B 


1000A 


7.0 


Gondhalekar et al. 1980 


C 


1000A 


4.9 






1400A 


5.0 






2200A 


2.5 


Henry et al. 1980 


D 


1400A 


5.3 






1180-1680A 


5.8 


Draine 1978 


E 


1000A 


6.7 






1400A 


9.8 






2200A 


3.7 



Notes to Table 2 

• Method A.- The FUV flux from bright stars B9 and earlier was inferred from its visual magnitude 
my and its color excess E(B — V). 

• Method B.- Similar to method A but using updated stellar atmospheres model and including only 
bright stars B5 and earlier. 

• Method C- Similar to method B but including more than 50,000 stars in the S2/68 Survey. 

• Method D.- Extrapolation from direct measurements of the spectrum in regions of 12° x 12° over 
about one-third of the sky. 

• Method E.- Values from the analytic expresion proposed by Draine (1978) for the "Standard UV 
background" . His formula is based on estimates of the typical FUV energy density at the solar circle 
and in the solar neighborhood (Habing (1968); Gondhalekar and Wilson (1975); Witt and Johnson 
(1973) and Jura (1974)) and on local measurements (Hayakawa, et al. (1969); Belyaev et al. (1971) 
and Henry et al. (1977)). 
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TABLE 3 

The Mean Energy Density at the Solar Circle 



Authors Method Wavelength (U) [10~ 17 ergcm~ 3 A -1 ] 



Habing 1968 A 


1000A 


3.9 - 6.7 




1400A 


3.4 - 5.6 




2200A 


1.4 - 2.3 


This paper B 


H2: [912-1100A] 


15.7 




FUV-H2: [1100-2070A] 


13.8 




FUV: [912-2070A] 


14.2 




1400A 


14.4 


This paper C 


H2: [912-1100A] 


15.9 - 11.4 




FUV-H2: [1100-2070A] 


12.1 - 9.8 




FUV: [912-2070A] 


12.7 - 10.0 




1400A 


12.5 - 10.3 



Notes to Table 3 

• Method A.- "Continuous method" (see §5.3). Contribution of 10 associations and 577 field stars 
between BO and B5. Left: direct star light; right: diffuse light included. Habing 1968, Table 6 and 
moderate extinction. 

• Method B.- As method A but using the PDMF (see §4.1). Diffuse light included (see §6, Appendix 
B). 

• Method C- From simulation. Left (right) values include (exclude) HII regions. Diffuse light included 
(see Sec. 6, Appendix B). The quoted values are the averages over 5 points 280 pc apart in a 
simulation of 57 Gyr. These mean values vary from point to point in one simulation and from one 
simulation to another (see text and Fig. 10). The variations of mean values excluding HII regions are 
much smaller. A point is inside an HII region about 2% of the simulation time. 
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TABLE 4 



The Typical Energy Density at the Solar Circle 



Authors 


Method 


Wavelength 


[7[10- 17 crgcm- 3 A- 1 l 

L O J 


Habing 1968 


A 


1000A 


2.5 < 4.0 < 14.0 






1400A 


2.5 < 5.0 < 15.0 






2200A 


1.5 < 3.0 < 13.0 


This paper 


median 


H2: [912-1100A] 


7.4 






FUV-H2- n 1 (10-2070 Al 


7.1 






FUV: [912-2070A] 


7.2 






1400A 


7.5 


This paper 


20% below 


H2: [912-1100A] 


4.6 






FUV-H2: [1100-2070A] 


5.0 






FUV: [912-2070A] 


5.0 






1400A 


5.4 


This paper 


20% above 


H2: [912-1100A] 


14.8 






FUV-H2: [1100-2070A] 


12.0 






FUV: [912-2070A] 


12.4 






1400A 


12.5 



Notes to Table 4 

• Method A.- Habing's 1968 recomended values for typical interstellar space near the solar circle (his 
Table 8). The three values given for each wavelength are the lower limit, the typical value, and the 
upper limit. 

• Median.- Median value of the field in a point in a simulation with t s i m — 57 Gyr. The first Gyr is 
discarded. See §4.2 and Figures 9 and 10. 

• 20% above.- The value below which the energy density exists 20% of the simulation time. 

• 20% over.- The value above which the energy density exists 20% of the simulation time. 



TABLE Bl 



Physical Properties of the Dust' 



Wavelength 




Kg 


9 


912-llOOA 


1.9 10~ 21 


8.95 10" 22 


0.8 


1100-2070A 


8.0 icr 22 


7.5 10" 22 


0.75 


1400A 


7.03 10" 22 


7.81 10" 22 


0.756 



a Units of k are cm 2 (H atom) 
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A. MASS DISTRIBUTION IN ASSOCIATIONS 



The IMF in equation (|) is the probability distribution for the stellar masses in an association. Each 
realization of an association will have fluctuations about this distribution, with the fluctuations being larger 
for smaller associations. In order to make a realization for a single generation of an association that is 
consistent with the IMF, we set the masses of stars with rm > 2Mq to 



(* = l,..,JV.) 



(Al) 



where 



■in 



i(Ph 



m h 



l/r 



(A2) 



and fio = m u . We have assumed that there are 5 generations in an association, so that the generation has 
A/"*fr,/5 stars. The mass m\ of the most massive star in each generation is chosen at random by setting 
x\ equal to a random number between and 1. The masses of the less massive stars are set assuming 
X2 = X3 = ... = 1/2. This procedure enables us to recover the correct IMF slope near m u when a 
distribution of associations is considered. Since an important fraction of £fuv an d Lh 2 is emitted by the 
most massive star in the association, this approach also accounts for the significant differences that can be 
present in similar (i.e. same M*h) associations. 



B. DIRECT AND DIFFUSE RADIATION 



In this appendix we describe the procedure used to calculate the scattered radiation, allowing for the 
possibility that the source and/or the observer are out of the galactic plane. We assume that the absorbing 
and dispersing medium follows the HI vertical distribution (eq. ^4|). The geometrical parameters of the 
problem are the height of the observer above the plane, z a , the height of the emitter, z e , and the projected 
separation between the two, x oe , as indicated in Figure Bl. 

To calculate the contribution dF sca .tt,i to the diffuse flux emitted at point 'e', scattered at point 'd', 
and observed at point 'o', we erect a coordinate system based on the line of sight 'oe' from the observer to 



the source. Let 'p' denote the projection of 'd' onto this axis, zu be the distance from 'p' to 'd' 
distance from the source to 'p'. We then have 



dK 



L * 



scatt,! 



4nr 2 



dfs, 



</>scatt(a,5) _ 7 



1 dzu vj d6 



and r ep the 



(Bl) 



where 



rf/scatt = k s n-a(zd) dr e 



nn(z d ) 



dr. 



(B2) 



is the fraction of the incident flux scattered per unit area (measured normal to the axis 'oe') at 'd'. In 
equation (B2), k s is the grain scattering cross section per H atom, r s is the scattering optical depth along 
the direct ray 'oe', and nn(zd) is the number density of H atoms at height = z p + w cos(9) cos(cj>) (see 
Fig. Bl) above the galactic plane. The term 4> SC att(ct, g) in equation (Bl) is the scattering phase function; 
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we adopt the Henyey & Greenstein (1941) phase function 

1 - g 2 

<t>Bca.tt(a, 9) = ,,-,. 2 o \3/2' ( B3 ) 

4 7r(l + g z - 2g cos a) A ' z 

where 

g=^-J dCl cos a<f) sca , tt (a,g)dQ (B4) 

is the average of the cosine of the scattering angle. The factors exp(— r e( j) and exp(— Tdo) in equation (Bl) 
are the extinctions along the paths 'ed' and 'do', respectively and include both absorption and scattering. 

For the 912 - llOOl and 1100 - 20701 bands, and at 1400A, the adopted values for n a , n s 
and g are given in Table Bl, and are derived from the data in Witt & Gordon (1999) assuming that 
n a = 5.34 x 10~ 22 (r/r y ) (1 - a) cm 2 (H atom)- 1 and k s = 5.34 x 10~ 22 (t/t v ) a cm 2 (H atom)" 1 . 

The total flux F sca tt,i that arrives at 'o' after being scattered once at 'd' is simply the integration of 
dF scat t,i over the whole volume. In practice, the integration is done over the revolution volume around the 
axis 'oe' that satisfies a < tt/2, because more than 96% of the diffuse flux comes from this volume. It is 
illustrative to compare -F sca tt,i to the direct flux 

F dil = -^e-^+ T °\ (B5) 



47rr 



oe 



where T a is the absorption optical depth along the ray 'oe'. Taking advantage of the fact that the scattering 
is strongly peaked in the forward direction (g > 0.75 — Table Bl), we normalize to the case of purely forward 
scattering (g = 1). In that case, the probability of a single scattering is r s exp(— r s ) and F scatti i = F d [ r T s . 
In the general case in which g < 1, we express the single scattering flux as: 

^scatt.l = F dil T s K. (B6) 

We have evaluated 1Z numerically in Figure B2. Note that it is close to unity, as expected from the fact 
that the scattering is strongly peaked in the forward direction. 

The procedure described above considers the diffuse radiation arriving at 'o' after only one scattering. 
To estimate the contribution of multiple scattering, we again consider the case of purely forward scattering. 
In this case, the probability that a photon will arrive at 'o' after n scatterings is r™ exp(— T s )/n\. We adopt 
the ansatz that the flux that arrives at 'o' after n scatterings when g < 1 is 

-Pscatt.n ~ -Pdir - S . - ■ (B7) 

n! 

The total flux F tot received at 'o' (direct plus scattered) is then 

F tot = F dir e^ = -^V ^ ~ Ts (1 " K) ■ (B8) 
4 7r r 

^ " 1 oe 

For purely forward scattering, 1Z = 1 and scattering has no effect on the observed flux, as expected. In the 
actual case, 1Z is close to unity and scattering has only a small effect. Figure B3 shows exp— r s (l — 1Z) as a 
function of the total optical depth for several different cases. 



The contribution of the diffuse radiation to the total radiation field obtained by this method can not 
be directly compared with the estimates obtained from observation because the latter only counts the 
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diffuse radiation coming from regions relatively far from sources. For example, the diffuse radiation that 
originates in lines of sight close to regions of recent star formation such as Gould's Belt are excluded. In §5, 
the median value of F sclitt /F tot from the simulations is about 0.2 whereas the estimates from observations 
are somewhat less than 0.1. However, if the scattered light that originates from lines of sight close to the 
sources are excluded from our calculations, our results are closer to the observations. We find that about 
1/3 of the single-scattered radiation has a scattering angle less than 5 deg; omitting such radiation from the 
"diffuse" category gives a diffuse fraction of about 13%, close to what is observed. 
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Fig. 1. — (a) The mean luminosity in the band [1100 — 2070 A] as a function of the stellar mass 
m. The triangles correspond to the calculated mean MS luminosities using stellar evolutionary 
tracks (see text) and the lines are the eight segment power law fit used to approximate the mass 
dependence of the luminosity in the simulations, (b) Same as (a) but for the H2 band [912, llOOA]. 
(c) The energy output in the FUV band [912 — 2070 A] for a mass distribution of stars following 
a Present Day Mass Function (PDMF) of MS stars. The thin lines correspond to the case when 
obscuration is taken into account, (d) Same as (c) but for the H2 band. 

Fig. 2. — (a) The mean main sequence ionizing photon production rate S49 (units 10 49 photons s _1 ) 
as a function of the initial stellar mass. The seven segment power law fit to these data is shown in 
broken line, (b) The contribution of the different mass ranges to the ionizing photon production 
rate assuming a Present Day Mass Function of MS stars. 

Fig. 3. — Various estimates of the IMF parameters in the plane [c^, — r] together with three curves 
of constant total ionizing photon production rate Sr,53 = Sr/(10 53 photons s _1 ). The letters beside 
the x symbols labels the estimates of Lequeux 1979 (L), Miller & Scalo 1979 (MS), Van Buren 1983 
(VB, b, c), and Scalo 1986 (Sc86). The 90% and 99% confidence contours correspond to Van Buren 
data. The cross near the label P indicates the adopted values of 4a and T in the present work. 

Fig. 4. — Evolution of the FUV luminosity for associations that ultimately form 1000, 100, and 10 
high-mass stars in five generations. The contribution by a single association to the energy density 
■ufuv {u-n = u/10~ 17 erg cm~ 3 A -1 ) at a point located 100 pc away, assuming that there is no 
extinction, is indicated in the right vertical axis. 

Fig. 5. — The projection on the plane of the positions of associations born in a period of 500 Myr 
(left panel) and 30 Myr (right panel). Large pentagons, large stars, small stars, and dots correspond 
respectively to associations with sizes characterized by A/"*,h stars with 1000 < A/"*^ < Af*h,u, 
100 < N*, h < 1000, 10 < N*,h < 100, and N* h ,l < N*,h < 10. The right panel reflects the OB 
associations, since the OB stars die after about 30 Myr. 

Fig. 6. — Expected increase of the mean energy density (E/Fuv)t sim i n the simulation time interval 
[0, t s im] as t S i m increases. Values normalized to the asymptotic mean value (C/fuv)oo- 

Fig. 7. — Time evolution of the FUV energy density at a point located at the center of simulation 
area and in the galactic plane. The bars in between the arrows indicate the periods of time where 
the point is inside an HII region; the labels give the characteristics of the association that produces 
the HII region: size A/"*,h, distance r, ionizing luminosity S49 and the mass mi of the most massive 
star. 

Fig. 8. — (a) Time evolution of the energy density in the [1100 — 2070A] band C/fuv-h 2 (continuous 
line) and in the [912 — llOOA] H2 band Un 2 (dotted line), (b) Ratio of the energy density in the 
[912 — llOOA] band, Uu 2 , to the energy density in the [1100 — 2070A] band, C/fuv-h 2 > as a function 
of the FUV energy density Up\jy. Each data point corresponds to the radiation field every 20 time 
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steps (i.e. 20 At a ~ 1 Myr) and the time period plotted is from t s - mi = 1000 to 4600 Myrs (see Fig 
7). (c) The fraction of U that comes to the observer as scattered light as function of U-n in the 
1100 — 2070A band. The times are the same as in Figure 8b. (d) The fraction of U in the FUV 
band that comes to the observer from the dominant source as function of U-n in the FUV band. 

Fig. 9. — (a) Distribution of the occurrence of C/h 2 values (dashed) and of C/fuv values (non- 
dashed). The arrows indicate the mean and median values of the distributions (see text §4.2). (b) 
Probability P(> U) that, at a given instant in the simulation, the FUV energy density is greater 
than U. The curves labeled 1 Gyr, 4 Gyr, and 56 Gyr, correspond respectively to the probabilities 
P(> U) calculated considering only the evolution of U in the time intervals At = 1 Gyrs, 4 Gyr, 
and 56 Gyr. The negative slope dashed line shows a power law fit to P(> U) for large values of 
U. The positive slope dashed line shows the corresponding power law for the typical time tjj (right 
vertical axis) between two events that raise the energy density to > U at a fixed point in space. 

Fig. 10. — Comparison of the evolution of mean energy density < U >t sim in a single simulation 
of time interval t s [ m and the expected relation (see Figure 6) of a large number of simulations. 
The corresponding evolution of the median value for a single simulation is also shown. Values 
normalized to the asymptotic mean value (J7fuv)oo- See text (§4.2). 

Fig. 11. — Evolution of the surface density of SN progenitors in the simulation. The continuous 
line corresponds to the average surface density over ^sim = 5 kpc 2 and the dotted line corresponds 
to the average within 500 pc of the central point. 

Fig. 12.— AU-n. run(^sim) (heavy line) and the total energy density C/fuv from the simulation 
(dotted line). The labels at the top give the properties of the runaway stars producing the most 
prominent FUV spikes (i.e. the stellar mass m, the velocity of ejection fejeo the distance of closest 
approach r m \ n to the point in question, and the generation in which the runaway was produced in 
the association). The labels at the bottom give the properties of the association where the runaway 
was born (i.e. the size of the association N*h and its distance to the observer rf orm ). 

Fig. 13. — Fig. Bl.- Geometrical parameters of the scattering problem. 

Fig. 14. — Fig. B2.- 1Z as a function of the distance r oe for z e and z D at the mid plane (continuous 
line), z = z e = 100 pc (long dash line), and z Q = 100 pc and z e = —100 pc (short dash line). 

Fig. 15. — Fig. B3.- Ratio of total flux to the flux for pure absobtion exp — r s (l — 1Z) as a function 
of the total optical depth for the three cases in Fig. B2. 
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